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iiXl'ANSION  AND  IMPROVFiMDNT  OI- 
SKF  COMPUTI-.R  PROGRAM  SHARP, RTII 


1.  - INTRODIIOTION 


The  objective  of  this  effort  was  to  produce  a design  and 
analysis  tool  in  the  form  of  a computer  program,  for  the  study 
of  ball  and  cylindrical  rolling  element  bearing  performance  under 
conditions  of  elastohydrodynamic  (EHD)  lubrication. 

The  program  developed  was  designed  to  treat  the  complete 
mill  t ibear  ing-shaft  system.  It  has  been  given  the  name  SllABPRTlI , 
an  acronym  for  "Sha f t -Rea ring -System-Thermal  Analysis". 

SllARIiRTIl  evolved  from  an  analysis  which  could  treat  highly 
loaded,  moderate  speed  shaft  bearing  systems  such  as  those  found 
in  helicopter  power  transmissions.  It  is  now  a design  and  analysis 
tool  addressing  simulation  of  moderately  loaded  but  ultra  high  speed 
gas  turbine  main  shaft  systems  anticipated  in  design  of  future 
aircraft.  The  latter  are  expected  to  operate  with  bearing  speeds 
greater  than  7>  million  DN.  fDN  is  the  product  of  the  bearing  bore 
diameter  expressed  in  millimeters,  and  the  shaft  speed  in  revolutions 
per  minute.)  The  difference  between  these  two  applications  noted  is 
the  degree  to  which  bearing  internal  friction  forces  affect  bearing 
performance . 

When  bearing  loads  are  high  and  speeds  moderate  the  bearing 
rolling  element -raceway  contact  traction  forces  are  large.  They 
have  the  potential  to  exceed  contact  inlet  drag,  as  well  as  cage 
friction  and  normal  forces  by  perhaps  an  order  of  magnitude.  At 
ultra  high  speeds,  under  moderate  loading,  the  relative  significance 
of  the  contact  friction  forces  decreases  with  respect  to  the  other 
friction  effects.  Rolling  element  and  cage  speeds  with  calculations 
based  on  epicyclic  or  raceway  control  assumptions  no  longer  reflect 
reasonable  estimates.  Thus,  the  dynamics  of  such  elements  must  be 
calculated  with  due  recognition  included  for  friction  and  inertia 
forces  which  act. 

The  increase  in  sophistication  required  to  incorporate  these 
calculations  in  the  description  of  rolling  element  and  cage  positions 
and  speeds  is  substantial.  Accurate  mathematical  models  of  the 
physical  phenomena  are  required.  Their  complexity  in  turn  requires 
state-of-the-art  numerical  mathematics  for  solving  the  resulting 
systems  of  highly  nonlinear  equations. 

To  an  extent  the  mathematical  friction  models  were  developed 
under  Air  Force  Contract  No.  F.X.361 5 -7  2 -C- 1467  and  Navy  MIPR  No. 
M62.376-.^-000007.  A cage  model  was  developed  under  NASA  Contract 
No.  NAS3-19739  and  installed  in  SllARFRTll.  Thus,  the  major  friction 
related  phenomena  had  been  modelled  and  could  be  used  successfully 
to  predict  bearing  friction  forces  and  heat  generation  rates 
present  in  heavily  loaded  bearings  in  which  cage  and  rolling 
element  speeds  could  be  estimated  with  conventional  methods. 
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The  lu'csent  effort  h;ul  to  address  and  incorporate  in 
SilABliRTH  refinements  in  modelling.  Additionally,  nonlinear 
etjuation  solution  techni(|ues  had  to  he  improvcii  to  treat 
the  am]ilified  influence  of  friction  forces  wliich  affect  perfor- 
mance in  moderately  loaded  high  speed  bearings.  Two  parallel 
approaches  were  taken. 


1.  The  first  was  to  investigate  and  improve  the  behavior 
of  the  m.athemat  i ca  1 models  in  SHABIiRTU.  The  equations 
which  they  create  were  to  be  as  linear  as  possible, 
while  remaining  consistent  with  physical  reality.  While 
making  these  investigations  methods  were  souglit  to 
improve  the  nonlinear  ecpiation  sol\'er,  snhV13,  present 
in  SIlABr.RTH. 

2.  The  second  effort  was  directed  at  investigating  several 
nonlinear  equation  solution  tcchni(]ues  which  had  a[ipearctl 
in  the  literature  and  were  believed  might  be  improved  to 
exceed  the  modified  Newton  Raphson- Fa  1 s i techniciues  of 
bOLVn. 

Details  of  tliese  two  efforts  shall  be  outlined  in  Sections  2 
and  7i. 

Under  this  contract  in  addition  to  enhancing  the  ability  of 
SHARF.RTH  to  produce  results,  its  actual  modelling  capabilities 
were  expanded  to  allow  bearing  ring  and  rolling  element  material 
properties  to  be  input,  individually.  Tliis  capability  will  allow 
proper  modelling  and  study  of  the  performance  of  bearings  having 
components  of  different  materials,  i.e.,  (^silicon  nitride).  The 
program  was  also  modified  to  permit  tlie  solution  of: 

1)  a single  cylindrical  roller  bearing  subject  to  radial 
and  moment  loading 

2)  a single  ball  bearing  subjected  to  axial  loading. 
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Bi'.ARiNC.  MOi)i:i.  MOD  ir  I CAT  IONS  anp  smAJTjON  sciiiiMi.  I mi*rovi;mi:nts 

Numerical  solution  for  rolling  element  bearing  quasi - 
dynamic  equilibrium  has  alwa^v  presented  difficulty.  A 
single  general  program  in  all  probability  will  not  be 
generated  to  solve  every  problem  to  any  desired  level  of 
accuracy.  However,  strides  have  been  made  to  obtain  such 
an  ideal  bearing  design  and  analysis  tool.  The  advent  of 
high  speed  digital  computers  has  enabled  continued  progress 
and  it  is  believed  that  progress  will  continue  as  long  as 
the  problem  is  pursued. 

The  computation  difficulty  resides  in  solving  large  sets  of 
nonlinear  equations  and  the  lack  of  .’..a  all  powerful  method 
for  that  purpose.  At  the  outset  of  this  portion  of  the 
effort,  it  was  believed  that  within  the  onstraints  of 
contract  time  and  money  the  modified  Regu  1 a- 1-a  1 s i , Newton- 
Raphson  iteration  techniques  of  bUl.VlS  were  tne  nest  avaiiaole. 
However,  recognizing  the  additional  computation  demands 
made  by  the  advanced  level  of  sophistication  in  this  work 
it  was  decided  to  attempt  an  upgrade  of  .S0I.V13.  The  improve- 
ments needed  could  be  obtained  by  recognition  of  the  following 
i deas : 

1)  The  solution  to  the  equation  set  would  be  easier  if 
the  equations  could  be  made 

a)  more  1 inear 

b)  load  variable-force  relationships  were  more  unique. 

2)  Achievement  of  a solution  would  be  more  jirobable  and 
faster  if  the  variable  values  were  confined  to  predeter- 
mined known  limits  within  which  they  must  lie. 

.3)  The  solution  would  be  faster  if  the  change  in  variable 
values  from  one  iteration  to  the  next  were  allowed  to 
be  only  within  certain  limits,  specified  uniquely 

for  each  variable  value.  Control  on  variable  value 
incrementation  is  ref._rred  to  as  damping. 

The  manner  in  which  these  ideas  have  been  pursued  are 
presented  below. 

2.1  Mathematical  Model  Changes. 

One  of  the  approaches  to  decreasing  the  difficulty  in  solv- 
ing nonlinear  equations  is  to  make  those  equations  less 
nonlinear.  This  approach  was  taken  in  the  revision  of 
both  the  rolling  element  cage  and  cage  ring  load  displace- 
ment models.  The  details  of  these  model  changes  are 
presented  in  the  Program  Users  Manual. 
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Another  substantial  change  in  the  mathematical  models 
was  the  change  in  the  concentrated  contact,  asperity 
friction  relationship.  This  calculation  was  changed 
from  a coulomb  model  to  a model  in  which  the  asperity 
friction  coefficient  is  a function  of  both^  contact 
sliding  and  rolling  speeds.  The  coulomb  model  requires 
the  contact  friction  force  to  be  the  product  of  a con- 
stant friction  coefficient  and  the  normal  load.  Thus, 
as  long  as  the  load  is  constant,  the  friction  force  can 
not  change  regardless  of  changes  in  rolling  element  speeds. 


The  Figure  1 demonstrates  the  difference  in  models 


Experience  with  SHABF.RTH  has  shown  that  asperity  traction, 
more  often  than  not,  is  responsible  for  larger  forces 
than  the  traction  resulting  from  lubricant  shear. 

2.2  Constraints  on  Variable  Values 

In  the  investigation  of  the  solution  of  a generally 
loaded  ball  bearing  test  problem  having  thirty-three 
degrees  of  freedom  (six  for  each  of  five  elements  and 
three  for  the  cage), it  was  observed  that  for  several 
sequential  iterations, a few  variables,  associated  with 
only  one  or  two  balls, were  the  most  influential.  Small 
changes  in  these  values  resulted  in  large  changes  in 
the  equilibrium  equation  residues.  As  the  influential 
variables  converged,  the  less  influential  variables 
drifted  to  impossible  values,  l/ltimately,  the  latter 
prevented  convergence  for  the  variable  set. 
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i The  variables  most  prone  to  the  above  behavior  were  the 

y and z components  of  ball  rotational  speed  in  lightly 
loaded  elements.  Both  of  these  components  take  on 
small  values  compared  to  the  x component  of  speed.  The 
nroper  sense  of  both  components  was  of  major 
importance  in  reaching  a converged  solution, 
i listabl  i shing  this  sense  must  be  done  prior  to  the  start 

of  the  solution  It  may  be  deduced  from  the  direction  of 
the  ou' er  and  inner  ring  rotational  speeds  and  the  rela- 
tive axial  position  of  the  inner  ring  with  respect  to  the 
outer  ring  at  the  particular  ball  location. 

The  logic  to  perform  the  above  function  was  introduced 
into  the  program  and  represents  a ref i nement  to  the  pre- 
vious version  of  the  code. 

2.7i  Solution  Damping 

The  Regula-Falsi , Newton-Raphson  scheme  in  SOI, VI. S, 
solves  a set  of  n nonlinear  equations.  The  right 
hand  sides  of  these  equations,  the  residues,  are  systema- 
tically reduced  and  approach  zero  at  solution.  The 
procedure  calculates  corrections  (AXin)  to  an  existing 
set  of  variables  (Xjisj)  according  to  I-q . 1. 

I 

^iN+1  " ^iN  ^^iN  i = l,...n 

where:  N indicates  the  current  itoratinn 

N+1  indicates  t!:c  next  iteration 
n is  the  tntal  TiumhoT-  q £ 'jrjkncv.TiS 


The  solution  is  said  to  be  converging  from  one  iteration 
to  the  next  if  each  set  of  corrections  has  the  affect 
of  reducing  the  equation  residues. 

In  the  solution  scheme  embodied  in  SOLVl,^,  Iic[ . fl)  is 
modified  as  follows: 


^iN+1  ^iN  ^iN  • ^^iN 


(21 


where:  is  a damping  factor  such  that 

For  each  iteration,  each  variable  has  an  independent 
damping  factor.  Four  steps  are  used  to  set  these  factors, 

1)  The  damping  factors  are  initialized 


f-M  = 1 
iN 


i = 1 


.n 


(31 


21  The  value  of  is  determined  using  Fq . (21, 

test  and  examined  according  to 

XMINI  .>X.j^^j>XMAXI  . (41 


If  either  test  is  true  is  halved,  hej . (2)  is  re- 

evaluated and  the  test  is  repeated.  This  process  is 
permitt^..d  to  occur  lOf)  times  for  each  variable. 

The  objective  of  this  scheme  is  to  restrict  the  value 
which  a variable  is  permitted,  to  that  sot  defined  by 

XMINl  .<  X.<Xf1AXI  . fS) 

After  step  2 the  individual  fj  values  may  be  different. 
XMIMj  and  XMAX I : are  preset  for  both  the  tliermal  and 
bearin}>  component  equilibrium  solutions. 

?>)  Tile  absolute  value  l^XiN/XiNl  is  calculated  for 
i=l,  n and  the  maximum  value  is  retained  and  compared 
to  the  maximum  value  \AX  jyi_  1 1 . If  the  value  at 

each  damping;  factor  is  modified  according  toP.q.  6. 


^^step3  ^^stepZ 


l<^^iN-l^^iN-  1 
I^XiN/XjN 


i = l 


(6) 


4)  This  step  was  taken  from  Wingo  (8).  At  this  point 
bq . (2)  is  evaluated  for  all  variables.  The  set 
of  equations  being  solved  is  then  evaluated  with 
XiN+1.  The  sum  of  squares  of  the  equation  residues 
is  calculated  and  the  test  in  Eq . (7)  is  performed, 

i=l  i=l 

If  the  test  in  Eq.  f7)  is  true  fj,  i=l,...n  is  halved  and 
step  4 is  repeated.  This  procedure  continues  fifteen 
times  or  until  the  test  fails. 

If  after  fifteen  halvings  of  fj,  the  sum  of  squares 
of  residues  is  not  reduced,  control  is  returned  to  tlie 
Newton  Raphson  scheme  without  any  correction  being 
applied  to  Xij,j.  The  step  size  used  to  evaluate  partial 
derivatives  is  altered  and  a new  set  of  i]s]  is  cal- 
culated. Steps  1)  through  4)  are  repeated.  As  many 
as  four  attempts  are  made  to  recalculate  ^Xj  in  a single 
Newton- Raphson  iteration  point.  If  the  test  in  Eq.  (7) 
has  not  been  passed,  the  scheme  terminates.  The 
variable  values  and  equation  residues  are  printed 
along  with  a warning  stating  "THIS  IS  THE  BEST  WE  CAN 
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DO  THIi  Rr.SULTS  MAY  RT,  USI-ABLi;".  Whether  or  not  tlie 
results  are  useable  depends  on  the  absolute  value  of 
the  magnitudes  of  the  residue  values.  Methods  to 
assess  useability  are  discussed  in  Ref.  11). 


Anothei  method  of  solution  damping  was  tested  in  place 
of  step  3)  above  and  although  it  worked  well  in  some 
instances,  it  actually  hindered  convergence  in  others. 
In  this  scheme  was  compared  to  Sj,  a vector 

whose  values  were  set  to  a constant  value  0.1,  or 
1.0  depending  upon  the  nature  of  Xj. 


AX/Xj  > Sj 


(«) 


If  for  a particular  variable,  i=l,  ...n  the  test  in 
liq . (8)  was  true  , for  that  variable  was  redefined 
such  that 


11’) 


This  scheme  thus  damps  the  corrections  to  the  variables 
individually.  It  was  found  that  when  two  or  more 
variables  were  closely  coupled,  convergence  was  pre- 
vented and  the  procedure  was  abandoned  in  favor  of 
step  3)  above. 
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Mi  llions  lOR  SOI.VIXC.  SYSTH''S  Or  NONMNIiAR  r.OIIATIONS 


The  sets  of  equations  generated  by  the  description  of  bear- 
ing-shaft interactions  in  SlIARf-RTH  are  highly  nonlinear  and  com- 
plex. The  increasing  complexity  in  detailed  representation  of 
the  physical  problem  has  necessitated  the  exploration  of  alter- 
nate solution  schemes. 

Tor  a set  of  N nonlinear  eepiations  with  \ variables,  Xj  , 

’>]';(  Xj)  = n i , i = 1 ,2 N'  (1) 

the  function 

= Z (2) 

i-l 

takes  on  the  minimum  value  11  = 0 at  all  points  satisfying  fq . f 1 ) . 

The  function  (1  defined  in  Tq . (2)  is  convenient  in  that  it  has 
the  same  form  as  the  least  squares  objective  function  in  regres- 
sion analysis  methods  for  generating  empirical  formulas  from  ex- 
perimental data.  Therefore,  some  of  the  w'ell  developed  solution 
techniques  in  regression  analysis  can  be  applied  to  solve  T.q . fl). 
•Additionally,  there  are  several  methods  available  for  calculating 
the  minimum  of  a function  of  many  variables.  These  methods  can 
also  be  used  to  solve  Fiq . (1)  expressed  in  the  form  of  Tq . (2). 

Two  methods  have  been  adapted  for  solving  T<i . (21,  viz., 
the  method  by  Tletcher  and  Powell  [2|  and  that  by  Rowell  [.1]. 

The  former  is  a steepest  descent  procedure.  It  uses  an  iterative 
[irocess  to  find  the  minimum  of  the  function  Ifi . The  basic  idea  is 
to  move  from  an  initial  point,  X,  along  tlie  vector  with  com.ponents 
[4) 

~ ' ' ■ 

whose  vaiues  change  continuously  as  the  path  is  followed.  The  meth- 
od described  by  Powell  is  a modified  version  of  the  (lauss  Newton 
technique.  It  involves  expanding  Tq . (1)  in  N Taylor  series  and 
uses  the  results  of  linear  least  squares  [4]  in  a succession  of 
i t e ra  t ions. 

In  general,  a steepest  descent  procedure  is  expected  to  con- 
verge fo:  poor  initial  values  but  requires  extensive  computation 
time.  The  f'.auss  Newton  technique,  however,  will  converge  rapidly 
for  good  starting  estimates.  In  the  ideal  situation,  use  of  the 


first  method  should  he  made  during  the  first  scv^eral  iterations. 

Then  a switcli  to  the  second  method  should  follow  to  achieve  rapid 
convergence.  The  solution  of  a set  of  nonlinear  eciuations,  however, 
remains  an  art  evolving  into  a science.  In  practical  situations, 
the  obtaining  of  a solution  is  more  often  than  not  related  to  the 
skill  and  past  experience  of  the  investigator  rather  than  guaranteed 
by  set  procedure  steps.  A clear  understanding  of  the  function  char- 
acteristics is  always  among  the  necessities  for  rapid  convergence 
no  matter  what  method  or  combination  of  methods  is  used.  Therefore, 
some  of  the  parameters  in  the  techniciues  discussed  have  to  be 
properly  constrained.  These  constraints  arc  dictated  by  the  fea- 
tures of  the  function  under  consideration  as  well  as  the  experience 
''ith  the  programs  on  the  part  of  the  investigator. 

RetcTrer  and  Powell  ^’ethod  .yy 

This  tTiethod  is  based  on  a theory  advanced  by  Davidson  [5|. 

T.xprcss  the  function  f in  the  quadratic  form 

> {A}.{X>  + i{X>.  [H]{X>  (-M 

where  ^ is  a constant  value  of  0 corresponding  to  the  assumed  or 
current  vector  {X>,  {A>  is  a vector  of  constants,  (X)  is  a solution 
vector  and  [II)  is  a matrix  whose  elements  arc  the  second  order 
partial  derivatives  of  ?)  with  respect  to  the  components  of  {X}. 

The  gradient  vector  {(',}  can  be  obtained  from  Ta) . (!S)  by 
differentiating  with  respect  to  the  components  of  {.X} , i.e. 

(C)  = (A>  + [111  {X>  (4) 

At  a minimum  value  of  0,  (X}  = {Xo}  and  {(ij  = 0.  Therefore, 

{A>  = -(ill  (Xo)  (5) 


or 


{(I)  = -[ll|{Xo>  * [1!1{X> 

Multiplying  both  sides  by  [II]  , one  obtains 

(Xo>  = {X)  - im'^  {r,>  (M 

In  this  method,  the  matrix  [11)"^  is  not  evaluated  directly 
but  approximated  in  the  iteration  process.  The  algorithm  proceeds 
as  fol lows : 


(11  Select  a starting  point  for  the  solution  vector  {X}. 
(2)  roiiipute  a tlirection  of  search  for  the  minimum, 


whe  re 


i^y 


[I- 


k)  _ 


[T] 


( k ) 


f> 


(k) 


(7) 


- 1 


K is  the  iteration  index  and  {S}  is  the  direction  vector.  It  is 
seen  from  F.q . (6)  that  {Sj  = ^Xo  -X>.  The  matrix  [ 1- ] can  initially 
be  represented  by  a diagonal  matrix  with  the  diagonal  elements  set 
equal  to  instead  of  the  unit  matrix  suggested  in  [1]. 

When  the  absolute  values  of  the  components  are  of  different 

orders  of  magnitude,  the  unit  matrix  can  seriously  impair  conver- 
gence and  in  actual  bearing  problems  lead  to  a convergence  stall. 
This  will  be  explained  in  the  next  step. 


chosen  hv 


A one  - il  i mens  i ona  1 
the  previous  step 


search  is  conducted  in  the  direction 
for  a minimum  utilizing  the  relation, 


{X} 

derived 


fk+n 


(k) 


(k) 


(8) 


from  Fci-s.  (6)  and  (7).  The  factor  ^ is  termed  the  step 
size.  When  close  to  the  solution  vector  {Xo>  , % is  unity  as  can 

from  Tqs.  (6,7).  The  ste|j).size  >\  is  calculated  so  that 
relative  minimum  at  A simple  algorithm  for  est- 

imating X can  be  found  in  f 4 f,  . Usual ly , the  relative 
exist  between  {Xj  and  ) 


be  seen 
0 has  a 


with  7S  in  I'q.  (81 


for 

mini mum  will 
expressed  by 


X = Minimum 


( 


-2H 


) 


fk) 


(9) 


However,  when  the  time  required  for  evaluating  the  function  is 
considerable,  it  is  advisable  to  use  the  value  of  7\  calculated 
according  to  Tq . (9)  as  long  as  the  function  value  of  decreases. 

The  algorithm  for  determining  ^ in  f4)  need  only  be  used  when  ^ 
diverges  or  becomes  very  small. 

It  is  noted  that  when  the  absolvite  valvies  of  the  comi'onents 
of  {X^  differ  by  orders  of  magnitude,  the  same  phenomenon  will 
will  exist  in  in  a hearing  problem.  If  a unit  matrix  is  assumed 

for  [I'l  in  I.q . (7),  then  {S}  = -{H},  and  X ‘■'■'•n  k>e  very  small  in 
b.().  (9)  for  a decreasing  function  0.  Thus,  the  improvements  in 
{X^  (k+1)  in  bc) . (8)  can  be  negligibly  small.  This  explains  the 
"convergence  stall"  mentioned  in  step  (2). 
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(4)  A convcrjicnce  check  hasci.1  on  a s[U'cifie(i  small  v.iliic 
ot  is  mailo.  If  convcrj’cncc  is  achicvctl,  the  procedure  is  ter 
n'.iiiated;  otherwise  a new  search  direction  is  chosen  per  steji  (2) 
with  [1|  calculated  as  follows: 


whe  re 


(k+1) 


(k) 


(A) 


(k) 


fH| 


( k ) 


( lOj 


;A| 


(k) 


IRl 


(k) 


{ Ax>n<Tt 

{a  /A  ni*  ^ I 

/ah  \Tkl  t r i;  I / . r A ) 


{ A X>  * = 

X 

1 

^,yk>l) 

in  which  t represents  the  transpose  of  a vector,  i.e.  a column 
vector  is  transposed  to  a row  vector.  A new  one  dimensional  search 
is  then  performed  in  the  new  direction.  The  process  is  repeated 
until  convergence  is  obtained. 


,3  2 Rowe  1 1 j^Iethod  .B 

The  basic  purpose  of  this  method  was  to  modify  the  (lauss 
Newton  technique  to  eliminate  the  necessity  of  derivative  calcula- 
tions. The  derivatives  are  approximated  by  finite  differences 
im[)licitly  without  requiring  substantially  more  function  evalu- 
ations. When  the  accuracy  in  calculating  derivatives  numerically 
is  not  dependable  and  the  number  of  function  evaluations  is  required 
to  be  minimal,  this  method  seems  to  be  very  attractive. 


b(is.  (1)  are  first  linearized  by  expanding  them  in  N Taylor 
^’Cries.  Only  linear  terms  are  retained, 


M 


'/'i  = <ki  + 21 


where 
va  1 lies 


' 4!’  * C5:)  , AX:  = Xj  - X3  X 

,'s  of  the  component.s  of  vector  fjji  . J 

Substituting  T.qs.  fll)  into  I'.q . (2)  and  then  setting 


(11) 

are  the  current 


1 ,2,  . . . ,N 


(12) 


1 1 


one  obtains  the  Gauss  Newton  formulation.  In  matrix  notation, 
this  can  be  written  as 


[M]  ^ [M]  {a 


113) 


whe  re 


[M]  is  the  Jacobian  matrix  with  elements  , t denotes  the 

transpose  of  a matrix  and  K is  the  iteration  index.  Note  if  ['1](k)t 
are  removed  from  both  sides  of  Eq . (13),  it  reduces  to  the  well 
known  Newton  Raphson  formulation. 

The  algorithm  of  the  method  proceeds  as  follows: 

(1)  A starting  point  for  the  solution  vector  {X>  is  selected 
and  a direction  vector  i^y  with  a component  of  unity  in 
each  coordinate  direction  is  assumed. 

(2)  The  Gauss  Newton  equations  are  set  up  and  solved  for 

(3)  The  vector  obtained  in  step  (2)  is  used  to  calculate 

a new  direction  vector,  which,  in  normalized  form,  can 

be  written  as 


(14) 


(4)  A one-dimensional  search  for  a relative  minimum  is  then 
conducted  in  the  new  direction  using  the  same  relationship 
represented  by  F.q.  (8).  The  procedures  described  by 
Powell  [ft]  are  used  to  find  the  step  size  X , which 

will  yield  a minimum  function  value  of  0 in  one  iteration. 

(5)  When  the  one-dimensional  minimum  has  been  found,  an 
overall  convergence  test  against  the  specified  small 
value  of  0 is  performed.  The  convergence  criterion 
suggested  in  [3]  is  that  when  both  (S>  and  X{S>  in  _ 

Eq.  (8)  have  acceptably  small  components,  then  iteration 
should  stop.  This  criterion,  however,  does  not  always 
guarantee  that  the  value  of  0 is  reasonably  small,  which 
is  necessary  for  convergence  to  be  established.  If  the 
convergence  criterion  is  satisfied,  the  procedure  is 


12 


terminated.  If  not, one  of  the  previous  direction  vectors  is  re- 
placed by  the  new  direction  vector.  The  vector  to  be  replaced  is 
the  one  with  index  corresponding  to  the  maximum  of  the  following 
quantities, 

where 

{P}  = [M]  ^ 


The  values  of  the  derivatives  in  the  f^’]  matrix  in  the  new 
direction  can  be  found  by  finite  differences  using  values  from 
the  completed  one-dimensional  search.  Thus,  the  fluass -Newton 
equations  are  now  updated  with  respect  to  the  new  direction  and 
are  solved  again  for  The  process  is  repeated  until  con- 

vergence is  achieved. 


It  is  noted  that  the  case  where  all  the  aforementioned  quanti- 
ties may  be  zeroes  is  not  considered  in  [-^].  When  it  happens, 
the  position  of  the  minimum  can  only  be  poorly  determined.  The 
situation  can  usually  be  corrected  by  changing  the  parameter  limit- 
ing the  step  size  which  is  to  be  supplied  to  the  computer  program 

[7]. 


3 . 3 Numerical  Px ample 


The  methods  discussed  have  been  put  in  a subroutine  subprogram 
named  BOVBER,  which  has  been  inserted  in  an  experimental  version 
of  SHABERTH.  Test  runs  have  been  made  on  a ball  bearing  problem. 
The  bearing  geometrical  dimensions  are  as  follows: 


(a) 

(b) 

(c) 

(d) 

(e) 


Pitch  diameter 
Diametral  clearance 
Ball  diameter 
Contact  angle 


190  mm 
0 mm 

23.0187  mm 
20  deg. 


Inner  and  outer  raceway  curvatures  0.51/0.52 


The  bearing  operating  conditions  are- 


Thrust  load  35580  Newtons 

Inner  ring  speed  13460  RPM 

output  contains  the  following  seven  unknowns: 


(a)  Ball  displacement  in  the  bearing  axial  direction. 

(b)  Ball  displacement  in  the  bearing  radial  direction. 

(c)  Ball  velocity  component  in  the  bearing  axial  direction, 

(d)  Ball  velocity  component  in  the  bearing  radial  direction. 

(e)  Ball  velocity  component  in  the  bearing  tangential  direction. 

(f)  Ball  orbital  velocity. 

(g)  Relative  displacement  between  ball  center  and  cage 
pocket  center. 


1 
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(•or  solving  the  set  of  seven  nonlinenr  eiju.'i  t i ons , the  methods 
discussed  i-'erfonned  very  poorly.  C 1 e t che  r- ;uui  - I’owe  1 1 ' s method 
reduced  the  value  of  (1  from  llfi.'l.S  to  ,'SO.  l in  7 iterations.  After 
that,  the  convergence  l-iecamc  very  slow.  The  comiiutation  was  then 
switched  to  I’owell's  methotl,  which  furtlier  reduced  0 to  18.7  in 
18  iterations,  and  then  experienced  a convergence  stall.  The 
problem  solver  suhroutine  .SOIA'i already  in  SHARI  inil  performed 
on  the  same  problem  with  excellent  results.  It  reduced  tlie  value 
of  (1  from  1185.8  to  2 . 4 8(i8X  1 0 ‘ in  merely  5 iterations. 


At  a fii'st  glance,  the  above  comparisons  clearly  indicate 
that  S0I.V1.5  is  superior  to  the  other  two  methods  exiilored.  Rut 
this  is  a misleading  generalization  from  particulars.  SOI, '15  has 
been  tested  extensively  on  bearing  problems  and  in  the  process 
lias  evolved  with  the  proper  constraints  for  rapid  convergence. 

Die  other  two  methods  have  not  undergone  the  same  refinements.  If 
the  three  methods  are  applied  to  a different  problem,  one  may  ex- 
pect different  results.  It  was  then  decide*'  to  solve  the  following 


Problem  wi 

i t h a 

set  of 

four  highly 

nonlinear  ecpiations. 

(Xi- 

X,)*- 

- (X,  - x^r 

* (Xi  - X2  - ^3  " 

^2  - 

X j s i n 

( 2 X 

.^1  X,  cos 

( 5 X,)  - 1 = 0 

f = 

^1  ^ 

7 

X“  . X 

3 . - 

3 -'4 

4 = 0 

4^4  = 

2X,  ^ 

3X3  ^ 4X^  - 

m = 0 

set  is 
1 

Assuming  starting  values  of  to  be 
Xi  = X,  = X3  = X^  = .5 

the  results  obtained  are  as  follows: 


A known  solution  for  this 


X,  = X, 


X = X = 
5 '4 


3.5.1  I- 1 etcher  and  Powell  Method 

In  this  method  the  initial  diagonal  matrix  for  |T|(sec  bq . (71) 
are  assumed  in  two  different  ways.  In  case  fa),  the  diagonal  ele- 
ments were  set  eipial  to 

in  the  first  four  iterations.  Succeeding  values  for  (1)  were 
calculated  using  lici.  flO).  Results  are  tabulated  as  shown  below: 
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VM 


(';isc  ( a ) 


No.  o i'  Iteration 

\'a  1 ue 

of 

1 

,4.0256 

x 

10^ 

2.4817 

X 

10^ 

3 

1.6101 

X 

10^ 

4 

2.041.4 

S 

1 . 52  45 

X 

10'^ 

(i 

1 . 8825 

X 

10'^ 

7 

4.1739 

X 

10'^ 

In  Case  fh) , the  same  assumption  for  the  diagonal  matrix  was 
used  only  in  the  first  iteration  and  a unit  matrix  used  in  the 
next  three  iterations  before  going  to  Eq . (10)  for  values  of  [E], 
Results  obtained  are  as  follovv's: 

(’  a s e ( h ) 

No . of  I terat  ion  Value  of  0 No.  of  Iteration  Value  of  f 


1 

4.0256  X 

10^ 

14 

2.3459 

X 

10--‘* 

2 

2.4817  X 

10'^ 

15 

1 .7255 

X 

1 0 - 

4 

2.0023  X 

10^ 

16 

1 . 1748 

X 

lo--'* 

4 

4.5718 

1 7 

7 . 8990 

X 

10-^ 

5 

1.4617 

18 

6.4554 

X 

10-' 

6 

2.6223  X 

10'^ 

19 

5.6346 

X 

10-' 

7 

8.9737  X 

10-2 

20 

4.6855 

X 

10-' 

8 

2.4277  X 

10"  2 

21 

2.5444 

X 

O 

1 

9 

1 . 4055  X 

10-2 

22 

7 . 4626 

X 

10-^ 

10 

8. 5634  X 

10-^ 

2 4 

2 .4264 

X 

10-^ 

1 1 

6.0581  X 

10--^ 

24 

1 . 5651 

X 

10  - ^ 

12 

4 . 5206  X 

10--^ 

25 

1 . 2606 

X 

10-^ 

13 

4.1359  X 

10--^ 

26 

8.7025 

X 

10-^ 

It  is  seen  that  the  cri'cct  of  the  assumed  values  for  [ I- | in 
the  first  \ iterations  has  great  influence  on  the  speed  of  conver- 
gence. Tt  is  believed  that  for  each  sjiecific  function,  there  is 
an  optimum  choice  for  [F],  F-rom  Fdis.  (7)  and  (8),  an  optimum 
choice  of  ^ should  also  speed  up  convergence. 

."S . 3 . 2 Powe  1 1 Method 

In  this  method,  the  input  data  include  a vector  w h i c h 

specifier  the  absolute  accuracy  limits  on  the  change  of  the  solu- 
tion vector  between  iterations  and  a parameter  liSCAI.F.  for 

limiting  the  step  size.  Different  values  were  assumed  for  {!  }■ 
and  FSC.M.F.  In  case  (a),  {F}  = 10‘4  and  FSCM.I  = 10-^.  Result 

are  tabulated  as  follows: 

C a s e (a  1 

No.  of  Iteration  Value  of  No.  of  Iteration  \’alue  of  ^ 


1 

3.0256 

X 

10^ 

14 

1 . 2850 

X 

. 7 

10  ‘■ 

2.6408 

X 

lo'^ 

15 

1 . 2308 

X 

7 

10  “ 

.3 

1 . 6054 

X 

10“^ 

16 

1 .2169 

X 

7 

10  ‘■ 

4 

4.9133 

X 

10^ 

17 

1 . 0051 

X 

7 

10-" 

.3 

1 .2008 

X 

10^- 

18 

7.6573 

X 

10-'^ 

6 

5.4462 

X 

10 

19 

6. 1980 

X 

10-'^ 

*7 

4.0010 

20 

4.6055 

X 

10-'^ 

8 

5.4152 

X 

10’  ^ 

21 

4.4769 

X 

10--^ 

9 

3.4274 

X 

10-2 

22 

2.3138 

X 

10-'^ 

10 

3.3390 

X 

10"2 

23 

1.6479 

X 

10-'^ 

1 1 

3.2167 

X 

10-2 

24 

3.2553 

X 

10-^ 

12 

1.4647 

X 

10-2 

2 5 

5. 1502 

X 

lo-'’ 

13 

1 .2851 

X 

10-2 

26 

2 .9703 

X 

10-2 

Total  no. 

of  function  evaluations  = 

1 32  . 

1 n C a s e 
results  were 

(h)  . {1}  = 

oht  a i ned : 

1 0 

" {X} 

and  IhSCATT.  = 10^. 

The  foil 

ow  i nji 

Case  ( h ) 

\'o.  of  Iteration 

\’a  1 ue 

of  0 

No . of  1 1 e rat i on 

Va  1 ue 

of  0 

1 

sC 
Ln 
r j 

X 

10^ 

1 2 

1.1717 

X 

1 o'  1 

*) 

2.5088 

X 

10^ 

13 

1.1184 

X 

10-' 

3 

2.4588 

X 

10^ 

14 

4 . 7794 

X 

in  “ 

4 

2.4550 

X 

10^ 

15 

8 .2500 

X 

10--^ 

5 

1.9918 

X 

lo'’ 

10 

5.  1900 

X 

10-^ 

h 

1 . 0901 

X 

10^ 

1 7 

1 . 8922 

X 

10'^ 

7 

1.2130 

X 

10^ 

18 

2 . 02  57 

X 

10'^ 

8 

4.9442 

X 

10^ 

1 9 

1.4810 

X 

10'^ 

9 

2 . 8403 

X 

10^ 

2 0 

1.7413 

X 

. 7 

10  ' 

10 

3. 1453 

21 

6.3507 

X 

lO"''' 

1 1 

9.3839 

X 

10'^ 

22 

1 . 0573 

X 

10'^ 

Total  no.  of  function  evaluations  = 9.3. 

In  case  fc),  {1.}  = in'‘^  {X}  and  FSCALH  = .3  X lO'^.  The 
results  obtained  are  shown  below: 


Case  (c) 


\o.  of  Iteration 

Va  1 lie 

of 

\o . of  1 1 e ra t i on 

Value 

of 

1 

3.0256 

X 

lo'^ 

14 

1 . 1683 

2 

2 . 5598 

X 

10^ 

1 5 

1 . 1 073 

3 

1.1546 

X 

lo'^ 

16 

1 . 1073 

4 

1 . 9929 

X 

10'^ 

17 

1 . 1073 

5 

7.7315 

X 

io2 

18 

9.2596 

X 

10'  ^ 

6 

7.2866 

X 

io2 

19 

9.2596 

X 

10'^ 

7 

5.7259 

X 

10^ 

20 

9.2596 

X 

10'^ 

8 

3.4214 

X 

10^ 

21 

9.2596 

X 

10'  ^ 

9 

3.4114 

X 

10^ 

22 

9.164  5 

X 

10'^ 

10 

3.0366 

X 

io2 

2 3 

8.1042 

X 

10'^ 

11 

2 . 0606 

X 

10 

24 

8.1042 

X 

10'  ^ 

12 

3.6928 

25 

8.1042 

X 

10'^ 

13 

1.2196 

26 

7.5514 

X 

10'' 

To^'al  no.  of  function  evaluations  = 151. 

It  is  seen  that  the  best  results  were  obtained  in  Case  (h)  . 
Increasing  the  value  of  F.SCALE  in  Case  (c)  resulted  in  failure  to 
converge.  These  results  clearly  indicate  that  there  are  optimum 
choices  for  these  parameters  for  any  specific  function  to  achieve 
rapid  convergence. 


I 


I 


I 

f 


3.3,3  S0I.V13 

1 he  solution  technique,  which  is  based  on  a conil>  i na  t i on  of 
the  Rcgula-Fa Isi  and  Newton- Raphson  methods,  was  used  to  solve 
the  same  problem.  The  results  calculated  follow: 


No.  of  Iteration  Value 

of  No.  of 

Iteration  Value 

_of 

1 

3. 02.36 

X 

lO'^ 

14 

8 . 6801 

x 

, 

1 0 

y 

2 . 9099 

X 

10^ 

1 5 

5.5226 

X 

10 

3 

1 . 8725 

X 

10^ 

16 

2 .9799 

X 

10 

4 

1 . 3936 

X 

10“^ 

1 7 

6.8156 

.3 

1.1845 

X 

10^ 

18 

4 . 8235 

. -> 

b 

7.0244 

X 

10- 

19 

2.2842 

X 

10  •" 

“7 

4.0705 

X 

10‘^ 

20 

2.9017 

X 

10'^ 

8 

2 . 3038 

X 

10'^ 

21 

6.2933 

X 

10-' 

9 

1.2547 

X 

10^ 

22 

1.1100 

X 

10"' 

10 

6.6686 

X 

10^^ 

2 3 

2.0917 

X 

lo"-'’ 

1 1 

3.2925 

X 

10^ 

24 

3.7563 

X 

10-^ 

12 

2.1160 

X 

10^ 

25 

5.7764 

X 

10'  ^ 

13 

1.4750 

X 

10^ 

From  the 

results  of 

the  second  problem. 

the  methods 

invest  i - 

gated  seem  to 

have  a definite  advantage 

over 

S0LV13,  e.g. 

Fletcher 

and  Powell.  It  is,  therefore , logical  to  expect  that  if  proper 
constraints  are  introduced  into  these  methods,  they  can  be  made 
into  very  effective  tools  for  solving  realistic  hearing  problems. 

The  proper  optimization  of  constraints,  however,  can  only  be  obtained 
through  experimentation  with  the  highly  nonlinear  equation  sets 
characteristic  of  bearing  analysis. 

As  the  detail  and  understanding  increase  in  computerized 
simulation  of  bearing  performance  current  equation  solvers  are 
taxed  to  the  limits  of  their  capabilities.  It  is  believed  that 
additional  experimentation  with  the  methods  described  will  yield 
the  more  powerful  tools  needed  to  meet  the  pressing  mathematical 
needs  of  state  of  the  art  design  and  analysis.  However,  current 
results  from  both  test  problems  highlight  tiie  development  needed 
for  the  alternate  schemes  evaluated,  Sf)LV13  is  therefore  the 
best  method  available.  It  is  thus  the  only  one  in  SIIARF.RTH 
at  this  time. 
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4.  - msciissioN  oi-  PunnuiMS 

4 . 1 I nt  rocluc  t ion 

Solutions  to  two  sample  problems  were  chosen  to  demonstrate 
the  use  of  SllABhRTII,  The  olijective  of  the  sample  problem  exe- 
cutions was  to  demonstrate  the  prime  ami  secondary  foptional) 
capabilities  of  SlIAbbRTll.  The  prime  capability,  to  calculate 
sha f t -bear i ng  system  performance  parameters  under  snecified 
conditions  of  load,  speed  and  temperature,  is  demonstrated  by 
the  analysis  of  the  high  speed  rolling  element  bearing  test  rig 
operated  by  the  Aero  Propulsion  Lalioratory,  Tubrication  Branch, 
at  the  Wr i ght - Pat t erson  Air  Force  Rase  . The 
study  of  the  rig  performance  also  demonstrates  the  optional 
calculation  of  bearing  operating  clearances  as  affected  by  system 
radial  thermal  gradients,  ring  rotational  sneeds,  rolling  clement- 
raceway loading  and  cold  shaft  and  housing  fits. 

The  sccomi  sample  problem  which  calculates  the  steady  state 
tliermal  performance  characteristics  of  the  CII-53  helicopter 
transmission  power  input  shaft,  demonstrates  the  SilABF.RTll  option 
to  calculate  system  steady  state  temperatures  as  a function  of 
bearing  frictional  heat  generation.  SlIARIiRTIl  can  also  calculate 
transient  thermal  performance  cliaractcr  i st  ics  and  has  b<  cn  used 
in  this  manner  to  predict  time  to  failure  of  a system,  subse- 
quent to  the  loss  of  lubrication.  Ref.  fPl. 

In  addition  to  demonstrating  the  optional  consideration 
of  physical  phenomena  which  affect  sliaft  bearing  system  per- 
formance, the  high  sj)eed  test  rig  problem  has  been  used  to 
generate  results  utilizing  NPASS=1,2  and  3 solution  levels  of 
SlIABIiRTil.  These  levels  consider  with  increasing  complexity, 
the  impact  of  friction  forces  upon  system  performance.  Computer 
solution  time  is  included  in  the  summary  of  results  to  demonstrate 
the  economics  of  each  solution  level. 

4 . 2 Discussion  of  the  Analysis  of  the  Wright-Pat terson  Aero 
Propulsion  Laboratory  Lubrication  Branch  High  Speed 
Rolling  Flcment  Rearing  Test  Rig 

Figure  2 is  a cross  section  of  the  modelled  test  rig. 

The  209  size  cylindrical  roller  bearing  labeled  L and  the  6220 
size  split  inner  ring  ball  bearing  labeled  A are  the  sliaft 
support  bearings  in  the  assembly.  A is  the  test  bearing  and  is 
shown  in  detail  in  Fig.  3.  All  input  data  to  SIIABFRTII,  required 
to  describe  bearing  A were  taken  from  Fig.  3 with  the  exception 
of  the  ball  and  raceway  asperity  slope  angles  wliicli  were  assumed 
to  be  two  degrees.  The  input  data  used  to  describe  bearing  L 
reflects  the  standard  SKF  Industries  Nil  209  configuration. 
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The  lubricant  type  ( MI  1-- L-  7 808G  ) , shaft  and  housing  fits, 
system  teinjie natures  fused  to  evaluate  bearing  operating 
clearances,  lubricant  v'iscosity  and  pressure  viscosity  coeffi- 
cient at  outer  and  inner  raceway-ball  contacts  and  for  the 
bulk  lubricant)  were  supplied  bv  the  l.ulir  i ca  t i on  Brandi  of  the 
Aero  Propulsion  haboratory.  The  estimates  for  the  percent  lubri- 
cant in  the  bearing  cavities,  the  film  replenishment  layer  thick- 
nesses and  the  asperity  friction  coefficients  were  made  by  SKP 
personnel  and  arc  consistent  with  the  recommendations  made  in 
the  SHARIiRTII  User's  Manual,  Ref.  (1). 

The  rig  in  Figure  2 appears  to  be  a three 
bearing  support  system.  However,  the  center  bearing,  denoted 
G,  is  loaded  axially  and  radially  through  its  housing  and  is 
the  means  by  which  the  shaft  is  loaded.  Thus,  for  SilARF.RTII 
input,  bearing  G is  replaced  with  the  radial  and  axial  shaft 
loading,  plus  a point  moment  estimated  to  arise  from  the  combined 
radial  and  axial  loading  on  G.  The  radial  and  axial  load  plus 
the  shaft  speed  were  provided  by  the  Lubrication  Branch. 

The  relatively  complex  shaft  geometry  is  well  witliin  the 
capabilities  of  the  flexible  shaft  analysis  in  SliABFRTIi.  Shaft 
internal  and  external  diameters  at  points  of  change  are  input 
and  arc  used  to  establish  the  shaft  flexural  cliaractci  istics . 


SilABF.RTlI  output  from  solution  level,  NPASS=2  is  presented 
in  Appendix  I.  The  output  should  be  self  explanatory.  However, 
the  SHABERTH  User's  Manual  provides  additional  detail  on  output 
input  definitions. 


SHABERTH  was  executed  with  similar  input  at  botli, levels  1 
and  3.  A comparison  of  results  and  computer  solution  time, 
provide  useful  insights  into  the  use'  of  SHABERTH  and  into  system 
performance  predictions. 


Solution  level  1,  solves  the  shaft-bearing  equilibrium 
equations  considering  rolling  element  raceway  elastic  contact 
forces  and  rolling  element  centrifugal  forces.  Ball  orbital 
and  rotational  speeds  are  estimated  using  methods  which  approach 
outer  raceway  control  theory  described  below. 


The  ball  speed  vector  pitch  angle  A is 
to  Eq.  r7.50),  Ref.  10  . ' 

tan  4 = — - - — 


calculated  according 


(15) 
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is  the  h.i  1 1 - on  t-e  r raceway  contact  anj^le 
= n/iln 

I)  is  the  hall  diarietcr 
dm  is  tlie  bearing;  pitcli  diameter. 

/4  corresponds  to  (IT-  tan  ^ f Ky/Wx ))  where  tan'^fWy/Wx) 
is  printed  in  SlIAhhRTIl  output. 

Given  ^ from  l':q . (15),  the  Wx  and  \Vy  components  of  ball 

rotational  speed  as  well  as  the  ball  orbital  speed  IVq  are 
calculated  assuninj;  no  relative  ball-raceway  slip  in  the  tangen- 
tial direction  at  the  center  of  eitlier  outer  or  inner  raceway 
contacts  . 

The  yaw  component  of  ball  rotational  speed  U'z,  is  estimated 

to  be 

IV z = -1.  x in  ■’  Ky  • Kq  (lb) 

where  is  the  hall  orbital  speed. 

Roller  orbital  and  rotational  speeds  arc  estimated  using 
epicyclic  conditions.  Upon  satisfying  this  shaft-bearing 
ecpiilibrium  condition  rolling  clement  orbital  and  rotational 
speeds  arc  recalculated,  taking  into  account  ring  displacements 
and  ball  contact  angle  variations.  A single  pass  is  then  made 
through  the  friction  and  lubrication  subroutines  wherein  clasto- 
hydrodynamic  (FillD)  film  thickness,  friction  forces  and  frictional 
heat  generation  rates  arc  calculated.  No  attempt  is  made  to 
arrive  at  that  set  of  rolling  element  positions  and  speeds  which 
will  place  each  element  and  the  cage  in  force  and  moment  eciuili- 
b r ium . 

Solution  level  NPASS=2  proceeds  as  docs  level  1,  through 
the  clastic,  system  solution.  The  relative  outer-inner  ring 
position  thus  obtained  is  held  fixed.  The  equilibrium  posi- 
tions and  rotational  and  orbital  speeds  for  all  rolling  elements 
and  the  cage  arc  determined  for  all  bearings,  one  bearing  at 
a time.  The  ecpiilibrium  state  requires  all  rolling  elements 
and  the  cage,  to  be  in  force  and  moment  equilibrium 
considering  all  rolling  element- raceway  and  rolling 
element-cage  normal  and  friction  forces.  Inertia  forces  including 
centrifugal  force  and  the  two  gyroscopic  moments  acting  on  cacli 
clement,  arc  also  included.  The  interaction  between  the  cage 
and  ring  on  which  it  is  piloted  is  considered.  Note  that,  with 


22 


71V 


I 


.r 


I 


I 

I 


the  frietion  forces  incliuletl,  the  sliaft  ajiplieci  loatiing  is  no 
longer  etiu  i 1 i Ii  ra  t eil  hy  the  hearing  react  ioti  forces. 

Tlie  \I’ASS=3  solution  proceeds  until  all  rolling  element 
and  cage  ecpjilihria  are  satisfieil.  Hearing  inner  ring  reaction 
forces  are  recpiired  to  equilibrate  the  shaft  api)lied  loading. 

At  this  level,  the  rolling  element  raceway  reaction  forces  which 
balance  tiie  applied  loading,  are  compriseil  of  friction  as  well 
as  normal  forces. 

Although  the  input  data  for  the  three  solutions  were  simi- 
lar they  were  not  identical.  In  both,  the  level  1 and  2 solu- 
tions, the  analysis  which  calculates  bearing  ope  rat  in/;  diametral 
c 1 earance , based  upon  cold  shaft  and  housing  fits,  thermal  gra- 
dients, ring  centrifugal  expansion  and  ring  loading,  was  employed. 
To  reduce  calculation  time  in  the  level  3 solution,  this  external, 
iterative  calculation  was  eliminated.  It  was  intended  to  specify 
the  operating  clearance  predicted  by  the  lower  level  solutions 
fin  both  cases,  -0.0228  mm  for  the  hall  bearing)  at  input.  The 
negative  number  indicates  a reduction  in  clearance.  Inadvertantly 
-0.00228  mm  was  actually  input.  The  only  sulistantial  impact 
this  had  on  results  was  to  increase  the  shaft  axial  displacement 
required  to  equilibrate  shaft  loading.  Table  1 shows  the  compari- 
son of  ball  bearing  results  from  the  three  solution  levels. 

The  roller  bearing  operated  very  close  to  epicyclic  conditions 
at  all  solution  levels  and  therefore,  that  data  is  not  tabulated. 

A fourth  set  of  results  is  presented  from  a second  NPASS=1 
solution  witli  the  estimates  of  the  ball  speed  vector  pitch 
angle  (as  calculated  by  F-q.  (15))  reduced  by  a factor  of  10. 

Also  Jiq.  (17),  rather  tlian  (16)  was  used  to  estimate  the  Wz 
component  of  rotational  speed 

W'7  = -0.75  Wv  Wy  • Wo 

^ |Wy  . Wo  I ^ ^ 

The  reasons  for  obtaining  this  second  NPASS=1  solution  are 
discussed  later. 

Before  the  results  of  Table  1 are  discussed,  it  's  impor- 
tant to  note  that  the  test  bearing  under  study  was  suLiected  to 
heavy  loads.  Additionally,  there  was  not  an  excessive  amount  of 
lubricant  in  either  bearing  nor  was  there  any  other  factor  which 
would  tend  to  induce  cage  slip.  Consequentially,  the  estimates 
of  rolling  element  orbital  speeds  used  in  the  level  1 solutions 
are  accurate  to  within  4%  of  the  level  3 solution.  As  a con- 
sequence of  the  low  pitch  angles,  and  the  absence  of  gross  sliding 
in  the  ball  raceway  contacts,  ball  orbital  speed  is  slightly 
higher  than  would  be  predicted  by  raceway  control.  This  effect 
accounts  for  the  4%  deviation  between  level  1 and  3 orbital 
speetl  predictions  notcii  above.  Should  significant  cage  slip  be 
predicted  in  the  level  2 solution,  level  1 results  would  differ 
substantially  from  level  2,  and  level  3 would  differ  from  level  2 
resul ts . 
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Tabic  1 points  out  that  for  the  test  rig  assembly,  the  level  1 
solution  provides  good  estimates  of  the  performance  of  the  system 
at  approximately  \A%  of  the  cost  of  the  level  2 solution  and  at 
approximately  1 , I'J  of  the  level  3 solution  cost. 

The  results  of  the  solution  level  .3  execution  represent  the 
accurate  solution  to  the  problem.  The  lower  level  solutions  are 
approximations.  For  estimates  of  hearing  deflections,  fatigue 
life  and  Idlli  film  characteristics,  all  solution  levels  produce 
similar  results.  Ihc  level  1 solution  using  raceway  control 
type  estimates  of  ball  speeds  is  at  variance  with  the  remaining 
runs  in  tlie  estimate  of  raceway  heat  generation  rates,  and  of 
course  tlic  estimate  of  component  speeds  themselves. 

ihc  purpose  of  the  solution  denoted  1*  was  to  demonstrate 
that  the  heat  rate  and  component  speed  estimates  at  level  1 can 
be  brought  into  line  with  the  level  2 and  3 solutions  in  a very 
simple  manner.  Only  two  program  statements  reqviire  change.  These 
changes  in  fact  were  not  made  to  fit  the  data  but  simply  reflect 
the  estimates  of  component  speeds  which  SHARF.RTIl  uses  to  begin 
the  level  2 and  level  3 solutions. 

The  component  speed  data  from  levels  2 and  3 indicates  that 
the  ball  speed  vector  becomes  substantially  more  parallel  with 
the  bearing  axis  than  would  be  predicted  by  raceway  control 
theory , 

The  low  pitch  angle  tendency  has  been  observed  often,  in 
many  different  ball  bearing  solutions  from  SHABERTll.  It  is 
postulated  that  insufficient  traction  forces  develop  in  the  ball 
race  contacts  to  resist  the  gyroscopic  moment  generated  by  ball 
rotation  about  an  axis  orthogonal  to  the  axis  of  ball  orbit. 

Mgz  = -JWoWy  (18) 

Typically,  ball  equilibrium  is  satisfied  with  relatively  small 
values  of  Wy.  This  observation  has  been  used  to  speed  solutions 
by  using  only  one  tenth  of  the  pitch  angle  calculated  with  Eq.  (16) 
as  the  initial  estimate. 

The  speed  vector  pitch  angle  predicted  by  SHABERTH  will 
increase  with  increasing  traction  coefficients.  This  has  been 
observed  using  a version  of  SHABERTH  containing  an  EHD  film 
tiiickness  model  developed  by  Lowenthal  et  al,  Ref.  (11)  and  a 
traction  model  developed  by  Allen,  Ref.  (12).  The  combination 
of  these  models  produce  higher  traction  coefficients,  than  do 
the  SKI-  models  in  SHABERTH.  The  identical  p?oblem  solved  with 
SHABERTH  defined  herein  and  the  alternate  version,  produced 
speed  vector  pitch  angles  of  practically  zero  and  greater  than 
would  be  predicted  by  outer  raceway  control,  respectively. 
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As  a matter  of  interest,  the  generality  of  the  SKf-  traction 
moilel  within  SIlAIUiRTII  provides  a very  simple  means  to  increase 
the  traction  coefficient  hy  simply  increasing  tlie  asperity 
traction  coefficient  specified  at  input.  To  date,  experience 
has  been  limited  to  use  of  a value  of  0.1  wrhich  is  consistent 
with  values  found  in  tiie  literature.  However,  Cheng  fl3)  re- 
cently has  suggested  that  values  up  to  approximately  0.5  may 
be  valid. 


The  prime  purpose  behind  the  sample  problem  executions 
was  to  demonstrate  a level  of  performance  capability.  The 
successful  execution  of  the  three  levels  of  solution  has  also 
shown  how  SHABIiRTII  should  be  used  most  economically  in  the  solu- 
tion of  a multi-bearing  system.  Namely,  if  a level  2 solution 
indicates  less  than  51  cage  slip,  the  level  1 solutions  will 
produce  good  estimates  of  bearing  deflections,  fatigue  life, 
and  l;lin  film  thickness.  The  level  2 solution  should  be  used 
to  make  an  accurate  assessment  of  bearing  heat  generation 
rates.  The  level  3 solution  should  be  used  only  when  cage  slip 
is  shown  to  be  significant  from  a level  2 solution.  When  level  3 
is  used  every  effort  should  be  made  to  simplify  the  problem. 
Temperature  maps  and  the  change  in  clearance  analysis  sliould 
not  be  executed  in  order  to  save  computing  costs. 

4 . 3 Discussion  of  CH-53  Steady  State  Modelling  Results 

The  output  from  the  Cll-53  test  problem  is  contained  in 
Appendix  2.  The  system  contains  one  216  size  cylindrical  roller 
and  a stack  of  four  6216  size  split  inner  ring  ball  bearings. 

The  bearing  heat  rates  were  calculated  executing  SHABF.RTH  at 
the  NPASS=1  level.  With  two  clearance  change  iterations  permitted 
for  each  thermal  iteration,  the  bearings  act  as  heat  sources 
which  supply  heat  to  the  temperature  nodes  representing  inner 
and  outer  rings,  rolling  elements  and  lubricant.  Additional 
heat  sources  include  a gear  mesh  and  a seal  which  arc  specified 
by  constant  heat  rate  values  at  input. 

All  initial  system  temperatures  were  set  at  an  estimated 
value  of  75°C.  Bearing  performance  is  then  determined  based 
upon  these  temperatures,  the  bearing  heat  rates  serve  as  input 
to  the  temperature  calculation  scheme  which  in  turn  produces  a 
new  set  of  system  temperatures.  These  new  temperatures  affect 
bearing  clearance  and  most  notably  lubricant  viscosity  in  the 
next  calculation  of  shaft  bearing  system  performance  parameters. 

Six  iterations,  calculations  of  bearing  system  performance.,  are 
required  before  the  equilibrium  condition  is  achieved.  Fquilibrium 


is  siitisfieil  when  a ciieck  of  all  system  tcmjierat  nres , comparing 
two  sequential  iterations  reflect  less  tlian  a 1”C  cliange. 

The  CII-53  noilal  map  and  the  ecpi  i 1 i h r i urn  temperatures  are 
presented  in  Fig.  The  numbers  which  appear  following  a slash, 

e.g.  /106  for  node  (7)  indicate  actual  measured  temperatures  obtained 
f ror  UeT . f 1 4 ) . 

Although  it  is  not  included  in  Appendix  2,  shaft-bearing 
system  output  data  may  optionally  be  printed  after  each 
calculation  of  bearing  frictional  heat  rates.  This  data  allows 
the  user  to  observe  system  performance  as  a function  of  tempera- 
ture. The  option  also  serves  as  a precaution  against  obtaining 
no  output  data  should  the  run  fail  in  a late  stage  of  execution. 

Since  the  program  input  is  printed,  the  thermal  data  input 
serves  as  a good  example,  to  be  reviewed  along  with  the  SHABFRTK 
User's  Manual,  when  preparing  thermal  data  for  another  system. 
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CONC!, US  IONS  AND  RliCOMMliNDAT  IONS 


S . - 


5 . 1 Cone  1 us i ons 


In  1972  an  effort  was  bcjjun  by  SKI-  Industries,  Inc.  under 
the  sponsorship  of  the  Air  force  and  Navy  to  develop  an  analytical 
tool  for  the  analysis  of  ball  bearings  in  gas  turbine  engines. 

The  prine  objective  of  this  effort  was  to  include  as  much  of 
the  knowledge  on  fill)  lubrication  as  possible,  in  a fashion 
easily  useable  by  a bearing  system  designer. 

With  the  completion  of  this  work  in  tlie  form  of  SKf 
Computer  Program  SHABfRTII,  not  only  the  ball  bearing  but  the 
cylindrical  roller  bearing  can  be  analyzed  with  state-of-the- 
art  lubrication  and  friction  models.  Additionally,  the  complete 
shaft  bearing  system  can  be  treated  with  the  complex  interaction 
of  shaft  and  bearing  displacements  properly  taken  into  account. 

The  effects  of  bearing  friction  with  full  or  partial  Hill) 
films  on  rolling  element  and  cage  dynamics  is  treated.  The 
effects  of  fill)  film  thickness  and  bearing  surface  topography 
on  bearing  fatigue  life  have  been  modelled  and  applied  to  both 
ball  and  cylindrical  roller  bearings.  The  bearing  cage  is 
modelled  in  a comprehensive  manner.  Rolling  element-cage  loading 
is  determined  as  a function  of  rolling  element-raceway  and  cage- 
ring friction  forces.  The  properties  of  the  shaft,  bearing  and 
housing  materials  may  be  input  thus  permitting  tlie  study  of  non 
steel  components. 


5.2  Recommendations 


fvaluation  of  the  physical  phenomena  considered  in  SlIARIiRTIl 
require  the  solutions  for  large  systems  of  nonlinear  equations. 

As  noted  in  this  report,  a reliable  generally  applicable,  tech- 
nique for  this  purpose  does  not  exist.  The  pursuit  of  such  a 
tool  should  be  encouraged  at  every  opportunity. 

SIIABI.RTII  in  its  current  form  offers  a unique  vehicle  for 
extending  the  capabilities  of  bearing  system  design.  The  follow- 
ing opj)ortuni  t ics  for  a cost-effective  return  on  additional 
design  and  analysis  tool  development  are: 
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n The  inclusion  of  tapered  and  spherical  roller  hearing 
modu 1 es . 

N’cecl:  Helicopter  transmissions  and  geared  turbo  fan  engines. 

2)  Ilxpansion  to  a multisliaft  system  with  shaft  interactions. 

Need : To  compute  inner  shaft  bearing  performance. 

31  A model  for  lubricant  distribution  witliin  a bearing 
which  would  treat  lubricant  delivery  methods  fjet  or 
through  race)  and  address  the  cage  as  an  active  rather 
than  passive  element. 

Need: U 1 tra  high  speed  bearing  operation  where  lubricant 
delivery  efficiency  is  a paramount  factor  in  obtaining 
the  required  operating  speed. 

4)  Rased  on  the  premise  that  bearing  instabilities,  such 
as  those  which  occur  in  sparsley  lubricated  gyro  bear- 
ings, result  from  changes  in  the  frictional  character- 
istics within  the  bearing,  SIlARhRTII  should  be  used  to 
evaluate  changes  in  bearing  internal  forces  which 
lead  to  instabilities. 

Need:  Gyro  bearing  instability  persists  as  a problem  with 
minimal  understanding  of  cause  and  effect  relationships. 

51  Add  a flexible  ring  analysis  to  the  cylindrical 
roller  module. 

Need: To  correctly  model  the  influence  of  out  of  round 
hearing  preload  when  this  approach  is  used  to  reduce 
cage  slip. 

61  The  analysis  should  be  expanded  to  permit  a flexible 
housing  represented  by  discrete  spring  constants  at 
each  bearing  location. 

Need : To  accurately  obtain  solutions  in  those  cages 
where  the  housing  stiffness  is  comparable  to  or  less 
than  that  of  the  shaft. 


7)  Create  rotor  - bearing  stability  analysis  module. 

Need : Tbe  nonlinear  displacement  load  response  of  bear- 
ings is  poorly  represented  in  current  stability  analysis 
programs.  The  detail  present  in  SIIABFRTII  would  create 
powerful  insight  in  the  general  stability  and  design  of 
rotor  bearing  systems  in  engine,  transmission  and 
machine  tool  applications. 

8)  Parametric  analysis  - program  exercise. 

Need : The  utility  and  optimum  return  on  sponsor  invest- 
ment in  as  complex  a program  as  SMABERTII  is  to  be  found 
by  two  procedures: 

a)  Specific  problem  simulation  - e.g.  steady  state  or 
time  transient  thermal  analysis  of  WPAFB  high  speed 
bearing  test  rig 

b)  Generation  of  an  enhanced  manual  [curves  and  guide- 
lines) - to  demonstrate  the  effects  of  various  input 
parameters  upon  system  performance.  This  will  reduce 
the  amount  of  program  usage  experimentation  which  each 
user  will  now  have  to  undertake. 

The  "pump-priming"  activity  of  b)  will  effect  wider 
recognition  for  SIIARFRTH  as  an  essential  tool  in 
the  state-of-the-art  bearing  system  creation  and 
evaluation. 

As  SIIABCRTH  is  used  throughout  the  industry,  additional 
ideas  for  enhancement  should  come  forth. 

SIIABFRTII  represents  a state-of-the-art  tool  for  the 
analysis  nf  hall  and  cylindrical  roller  hearing  shaft  systems. 

Its  use  should  be  encouraged  within  the  industry  to  upgrade  bearing 
system  design  and  to  provide  a degree  of  commonality  in  bearing 
calculations  among  mechanical  system  suppliers. 
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FIGURE  3.  6220  SIZE  SPLIT  INNER  RING  TEST  BEARING 
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WPAFR  LUBRICATION  BRANCH  ROLLING  ELEMENT 
BEARING  TEST  MACFIINE  ASSEMBLY  I 

SHABERTH  OUTPUT 
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-P1F3  LU39IC4TI0N  BRANCH  • R.E.8.  TEST  MACHINE  ASSEHScr  I • SOLUTION  LEUEL 
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